program VMD

  use fft_wrap, only   : fft_init
  use mpi_wrap
  use efield
  use evolve
  use collisions

  use dimensions !, only : jx, jz, jv, jDir

  use parameters, only : restart, nt, dt, iDiag, iData, vPhase, runname
  use parameters, only : iTransAdvect, iEraseNonlin, iLinVlasov, iLinHydro, iFourier

  implicit none

  !-- major arrays ---------------------

  complex, allocatable  :: dfu(:, :, :, :, :)   ! distribution function
  complex, allocatable  :: uxMat(:, :, :, :)    ! velocity, x-component
  complex, allocatable  :: uyMat(:, :, :, :)    ! velocity, y-component, unused
  complex, allocatable  :: uzMat(:, :, :, :)    ! velocity, z-component, unused
  complex, allocatable  :: Ex2D(:, :, :)        ! electric field, x-component 
  complex, allocatable  :: Ey2D(:, :, :)        ! electric field, y-component 
  complex, allocatable  :: Ez2D(:, :, :)        ! electric field, z-component

  real,    allocatable  :: filter(:,:,:)        ! conical filter
 
  real,    allocatable  :: f0vec(:)             ! background distribution, gaussian
  complex, allocatable  :: fInit(:,:,:)         ! initial df  (1D equilibrium)


  !-- variables  ----------------------- 

  real       :: t                               ! time
  real       :: courant
  real       :: Courant_in_x                    ! advection Courant number 
  real       :: Courant_in_y                    ! advection Courant number 
  real       :: Courant_in_z                    ! advection Courant number 

  integer    :: it                              ! iteration count
  integer    :: iargc                           ! argument count


!-------------------------------------------------------------------------------------

  allocate ( dfu(jv, jz, jy, jx, jDir) )
  allocate ( uxMat(jz, jy, jx, jDir) )
  allocate ( uyMat(jz, jy, jx, jDir) )
  allocate ( uzMat(jz, jy, jx, jDir) )
  allocate ( Ex2D(jz, jy, jx) )
  allocate ( Ey2D(jz, jy, jx) )
  allocate ( Ez2D(jz, jy, jx) )

  allocate ( filter(jz,jy,jx) )

  allocate ( f0vec(jv) )
  allocate ( fInit(jv,jz,jDir) )

 
      !--  reading the run_name from the command line

      if(iargc().eq.0) then
         call getarg(0,runname)
         write(*,*)      'Usage: ',trim(runname),' <run name>'
         call mpi_abort(MPI_COMM_WORLD, -1, mpi_err)
         stop
      end if
      call getarg(1,runname)


      !-- set parameters, coordinate space, and the filter --

      call mpi_wrap_init()

      call read_input()

      call fft_init(dfu)

      call set_space_vectors()
      call set_filter(filter)

      filter(:,:,:) = 1.0      ! uncomment for comparibility with BGK code

      call evolve_init()
      
      call efield_init()



      !-- restore from restart file or create initial conditions --

      call set_ic( dfu, uxMat, uyMat, f0vec, fInit, it)

      t   = it*dt


      call efield_evolve_E0(t)

      call efield_compute_E(dfu, Ex2D, Ey2D, Ez2D, filter)


      !-- write initial diagnostics and binaries --

      if (.not. restart) then

         Courant_in_z = vPhase*dt/dz

         call diagnostics( dfu, Ez2D, Ex2D, Ey2D, f0vec, t)

         !call evolve_Uz( dfu, f0vec, uzMat)
         call write_bin(dfu, uxMat, uyMat, it, it/iData)

         write(msg, "('Output ' i4.4)") it/iData
         call mpi_wrap_msg

      end if

      write(msg,*) char(10), ' Begin temporal evolution', char(10)
      call mpi_wrap_msg


!-------------------------------------------------------------------------------------
! BEGIN EVOLUTION
!-------------------------------------------------------------------------------------

call evolve_F_collisions(dfu, dt/2)

do while (it .le. nt) 

   call dealias_5D(dfu)

   call evolve_F_long_advect(dfu, .true.)

   call efield_evolve_E0(t)
   call efield_compute_E(dfu, Ex2D, Ey2D, Ez2D, filter)
   call evolve_F_accel(dfu, Ez2D, f0vec)

   call evolve_F_long_advect(dfu, .true.)

   t  = t + dt
   it = it + 1
   
   !-- diagnostics --

   if ( (it/iDiag)*iDiag.eq.it ) then      
   
      call evolve_F_collisions(dfu, dt/2)

      call efield_evolve_E0(t)
      call efield_compute_E(dfu, Ex2D, Ey2D, Ez2D, filter)

      call diagnostics( dfu, Ez2D, Ex2D, Ey2D, f0vec, t)

      !write(msg, "(' Diagnostcs, t =' e12.6)") t
      !call mpi_wrap_msg
     
      
      if ( (it/iData)*iData.eq.it ) then

         !call evolve_Uz(dfu, f0vec, uzMat)
         call write_bin(dfu, uxMat, uyMat, it, it/iData)
         write(msg, "(' Output ' i4.4)") it/iData
         call mpi_wrap_msg

      endif

      call evolve_F_collisions(dfu, dt/2)

   else
      call evolve_F_collisions(dfu, dt)
   endif

end do      ! END TEMPORAL LOOP

    write(msg,*) char(10), ' End temporal evolution', char(10)
    call mpi_wrap_msg

    call mpi_wrap_exit()
    stop;

end


!-------------------------------------------------------------------------------------
